Packet 7 dimensionality reduction

library(tidyverse)
library(langcog)
library(psych)
library(readxl)
library(cowplot)
library(lme4)
library(lmerTest)
library(kableExtra)
library(lubridate)
library(ggdendro)
theme_set(theme_bw())

heatmap_fun <- function(efa, factor_names = NA){
  
  # get factor names
  if(is.na(factor_names)){
    factor_names <- paste("Factor", 1:efa$factors)
  }
  
  # put factors in a standard order when applicable
  body_factors <- factor_names[grepl("BODY", factor_names)]
  
  leftovers <- factor_names[!factor_names %in% body_factors]
  heart_factors <- leftovers[grepl("HEART", leftovers)]
  
  leftovers <- leftovers[!leftovers %in% heart_factors]
  mind_factors <- leftovers[grepl("MIND", leftovers)]
  
  other_factors <- leftovers[!leftovers %in% mind_factors]
  
  factor_levels <- c(body_factors, heart_factors, mind_factors, other_factors)
  
  # get factor loadings
  loadings <- efa$loadings[] %>%
    data.frame() %>%
    rownames_to_column("capacity") %>%
    gather(factor, loading, -capacity) %>%
    mutate(factor = as.character(factor(factor, labels = factor_names)),
           factor = factor(factor, levels = factor_levels))
  
  # get fa.sort() order
  order <- loadings %>%
    group_by(capacity) %>%
    top_n(1, abs(loading)) %>%
    ungroup() %>%
    arrange(desc(factor), abs(loading)) %>%
    mutate(order = 1:length(levels(factor(loadings$capacity)))) %>%
    select(capacity, order)
  
  # get percent shared variance explained
  shared_var <- efa$Vaccounted %>%
    data.frame() %>%
    rownames_to_column("stat") %>%
    filter(stat == "Proportion Explained") %>%
    select(-stat) %>%
    gather(factor, var) %>%
    mutate(factor = as.character(factor(factor, labels = factor_names)),
           factor = factor(factor, levels = factor_levels)) %>%
    mutate(var_shared = paste0(factor, "\n", round(var, 2)*100, "% shared var.,"))
  
  # get percent total variance explained
  total_var <- efa$Vaccounted %>%
    data.frame() %>%
    rownames_to_column("stat") %>%
    filter(stat == "Proportion Var") %>%
    select(-stat) %>%
    gather(factor, var) %>%
    mutate(factor = as.character(factor(factor, labels = factor_names)),
           factor = factor(factor, levels = factor_levels)) %>%
    mutate(var_total = paste0(round(var, 2)*100, "% total var."))
  
  # make plot
  plot <- ggplot(loadings %>% 
                   left_join(order) %>%
                   left_join(shared_var %>% select(-var)) %>%
                   left_join(total_var %>% select(-var)) %>%
                   mutate(capacity = gsub("_", " ", capacity),
                          factor = factor(factor, levels = factor_levels),
                          xlab = paste(var_shared, var_total, sep = "\n")),
                 aes(x = reorder(xlab, as.numeric(factor)), 
                     y = reorder(capacity, order), 
                     fill = loading, 
                     label = format(round(loading, 2), nsmall = 2))) +
    geom_tile(color = "black") +
    geom_text(size = 3) +
    scale_fill_distiller(limits = c(-1, 1), 
                         palette = "RdYlBu",
                         guide = guide_colorbar(barheight = 10)) +
    theme_minimal() +
    scale_x_discrete(position = "top") +
    theme(axis.title = element_blank())
  
  return(plot)
  
}
s_moments <- function(p) {p*(p+1)/2}
param_est <- function(p, k) {p*k + p - (k*(k-1)/2)}
check_ok <- function(p, k) {
  a <- (p-k)^2
  b <- p+k
  return(ifelse(a>b, TRUE, FALSE))
}
max_ok <- function(p) {
  df_check <- data.frame()
  for(i in 1:p){
    df_check[i,"check"] <- check_ok(p,i)
  }
  max <- df_check %>% filter(check) %>% nrow()
  return(max)
}
reten_fun <- function(df, rot_type = c("oblimin", "varimax", "none")){
  
  # figure out max number of factors to retain
  n_var <- length(names(df))
  max_k <- max_ok(n_var)
  
  # run efa with max factors, unrotated
  fa_unrot <- fa(df, nfactors = max_k, rotate = "none", 
                 scores = "tenBerge", impute = "median")
  eigen <- fa_unrot$Vaccounted %>%
    data.frame() %>%
    rownames_to_column("param") %>%
    gather(factor, value, -param) %>%
    spread(param, value) %>%
    filter(`SS loadings` > 1, `Proportion Explained` > 0.05)
  retain_k <- nrow(eigen)
  
  fa_rot <- fa(df, nfactors = retain_k, rotate = rot_type,
               scores = "tenBerge", impute = "median")
  
  loadings <- fa_rot$loadings[] %>%
    data.frame() %>%
    rownames_to_column("capacity") %>%
    gather(factor, loading, -capacity) %>%
    group_by(capacity) %>%
    top_n(1, abs(loading)) %>%
    ungroup() %>%
    count(factor)
  retain_k_final <- nrow(loadings)
  
  return(retain_k_final)
}
source("./scripts/p7_data_prep.R")

Joining, by = c("p7_entr", "p7_2day", "p7_ver", "p7_batc", "p7_resample", "p7_ctry", "p7_subj", "p7_file", "p7_recr", "p7_wher", "p7_whoc", "p7_abs_child.exp", "p7_abs_poetic", "p7_abs_tv.real", "p7_abs_see.image", "p7_abs_mind.world", "p7_abs_clouds", "p7_abs_vivid.dreams", "p7_abs_mystic.exp", "p7_abs_step.outside", "p7_abs_textures", "p7_abs_too.real", "p7_abs_music.attn", "p7_abs_heavy.body", "p7_abs_sense.presc", "p7_abs_fire", "p7_abs_nature.art", "p7_abs_colors", "p7_abs_thght.wander", "p7_abs_vivid.past", "p7_abs_makes.sense", "p7_abs_become.chctr", "p7_abs_visual.thghts", "p7_abs_small.things", "p7_abs_music.lift", "p7_abs_noise.music", "p7_abs_scented.mem", "p7_abs_visual.music", "p7_abs_before.said", "p7_abs_physical.mem", "p7_abs_voice.sound", "p7_abs_not.physical", "p7_abs_thgts.image", "p7_abs_odor.to.color", "p7_abs_sunset", "p7_abs_total", "p7_abs_check", "p7_dse_god.prescn", "p7_dse_conect.life", "p7_dse_no.daily.conc", "p7_dse_spi.strength", "p7_dse_spirt.comfort", "p7_dse_inner.peace", "p7_dse_god.help", "p7_dse_guided.daily", "p7_dse_direct.love", "p7_dse_lov.thru.othr", "p7_dse_touch.by.beau", "p7_dse_blessings", "p7_dse_selfless.care", "p7_dse_accept.wrong", "p7_dse_total", "p7_dse_check", "p7_se_voice.out", "p7_se_voice.in", "p7_se_placed.thought", "p7_se_vision.out", "p7_se_image.in", "p7_se_touch", "p7_se_smell", "p7_se_taste", "p7_se_dream.sent", "p7_se_stand.beside", "p7_se_demon.in.room", "p7_se_spnat.presence", "p7_se_shaking.prayer", "p7_se_emotion.prayer", "p7_se_powrful.prayer", "p7_se_out.body.exp", "p7_se_body.control", "p7_se_slep.paralysis", "p7_se_god.thru.pain", "p7_se_god.illness", "p7_se_live.healing", "p7_se_own.healing", "p7_se_total", "p7_se_check", "p7_wob_set.mind_reverse", "p7_wob_find.ways_reverse", "p7_wob_own.hands_reverse", "p7_wob_future.on.me_reverse", "p7_wob_little.change", "p7_wob_helpless", "p7_wob_others.do", "p7_wob_beynd.control", "p7_wob_interfere", "p7_wob_little.contrl", "p7_wob_cant.solve", "p7_wob_pushed.around", "p7_wob_total", "p7_unev_voice.aloud", "p7_unev_phone.ring", "p7_unev_call.name", "p7_unev_music", "p7_unev_no.ones.vox", "p7_unev_shadows", "p7_unev_total", "p7_unev_check", "p7_exsen_esp.exists", "p7_exsen_esp.exp", "p7_exsen_psychic", "p7_exsen_view.future", "p7_exsen_dream.true", "p7_exsen_dist.msg", "p7_exsen_send.msg", "p7_exsen_total", "p7_exsen_check", "p7_hthk_complex", "p7_hthk_responsblt", "p7_hthk_not.fun", "p7_hthk_lil.challeng", "p7_hthk_avoid.think", "p7_hthk_long.hrs", "p7_hthk_hrd.hav.to", "p7_hthk_smal.daily", "p7_hthk_lil.thought", "p7_hthk_way.to.top", "p7_hthk_new.soltions", "p7_hthk_not.exciting", "p7_hthk_puzzles", "p7_hthk_abstract", "p7_hthk_intel.task", "p7_hthk_mental.effrt", "p7_hthk_job.done", "p7_hthk_not.personal", "p7_hthk_total", "p7_por_thgs.hrt", "p7_por_thgs.hurt_a", "p7_por_thgs.hurt_b", "p7_por_thgs.hurt_c", "p7_por_wifi.thgs", "p7_por_job.wish", "p7_por_angr.cntrl", "p7_por_angr.cntrl_a", "p7_por_angr.cntrl_b", "p7_por_angr.cntrl_c", "p7_por_sprt.envy", "p7_por_read.thgs", "p7_por_read.thgs_a", "p7_por_read.thgs_b", "p7_por_read.thgs_c", "p7_por_stre.spoil", "p7_por_stre.spoil_a", "p7_por_stre.spoil_b", "p7_por_stre.spoil_c", "p7_por_conslt.unseen", "p7_por_mircl.prayer", "p7_por_pry.dead.back", "p7_por_spkn.curse", "p7_por_spkn.curse_a", "p7_por_spkn.curse_b", "p7_por_spkn.curse_c", "p7_por_curse.sick", "p7_por_sprt.put.thgs", "p7_por_fall.in.lov", "p7_por_fall.in.lov_a", "p7_por_fall.in.lov_b", "p7_por_fall.in.lov_c", "p7_por_thgs.heal", "p7_por_visualization", "p7_por_total", "p7_por_check", "p7_mm_ang_feel.hurt", "p7_mm_ang_thgs.hurt", "p7_mm_ang_sprt.hurt", "p7_mm_ang_physical", "p7_mm_ang_sickness", "p7_mm_car_fel.no.pr", "p7_mm_car_thk.no.pr", "p7_mm_car_sprt.help", "p7_mm_car_physical", "p7_mm_car_curing", "p7_mm_env_feel.hurt", "p7_mm_env_thgs.hurt", "p7_mm_env_sprt.hurt", "p7_mm_env_physical", "p7_mm_env_sickness", "p7_mm_thnk.feel.hurt", "p7_mm_sprt.thgs.hurt", "p7_mm_total", "p7_mm_check", "p7_dem_sex", "p7_dem_age", "p7_dem_pocc", "p7_dem_major", "p7_dem_ethnicity", "p7_dem_rur.urb", "p7_dem_affrd.basics", "p7_dem_ses", "p7_dem_how.sprt.relg", "p7_dem_religion", "p7_dem_church", "p7_dem_holy.tung.gif")
NAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercion

d1 <- d0 %>%
  select(-contains("_dem_"), -ends_with("_cat"), -p7_ctry,
         -contains("total"), -contains("check")) %>%
  gather(question, response, -p7_subj) %>%
  # rescale everything to be in 0, 1
  mutate(scale = gsub("p7_", "", question),
         scale = gsub("_.*$", "", scale),
         response = case_when(
           scale %in% c("abs", "exsen") ~ response,
           scale == "dse" ~ response / 5,
           scale == "hthk" ~ (response + 2) / 4,
           scale %in% c("mm", "unev") ~ response / 3,
           scale == "por" ~ response / 2,
           scale == "se" ~ response / 4,
           scale == "wob" ~ (response + 3) / 6)) %>%
  # get rid of follow-up questions for porosity
  filter(!grepl("_.$", question)) %>%
  select(-scale) %>%
  spread(question, response) %>%
  column_to_rownames("p7_subj")

This is an exploration of dimensionality reduction on Packet 7 (last updated: 2019-05-19). KW thought of this in response to the general question of how (if at all) we can distinguish porosity “beliefs” from spiritual “experiences.”

The basic question here is whether we can distinguish beliefs from experience by tracking patterns of covariance across individual participants: What are the “clusters” of questions that tended to hang together at the participant level (i.e., if someone endorsed one question, what else did they endorse)? If porosity and spiritual experience are fully redundant with each other, we would not expect to see separate clusters of porosity vs. spiritual experience items.

Main take-away: If we retain more than 2-3 factors, we end up distinguishing between porosity and spiritual experience. Two standard methods of determining how many factors to retain (parallel analysis and minimizing BIC) both suggest retaining six or more factors.

Exploratory factor analysis (EFA)

Parallel analysis

fa.parallel(d1)

Parallel analysis suggests that the number of factors =  12  and the number of components =  11 

Parallel analysis suggests retaining 13 factors.

efa13 <- fa(d1, 13, rotate = "varimax")

heatmap_fun(efa13, factor_names = paste0("MR", 1:13))

the condition has length > 1 and only the first element will be usedJoining, by = "capacity"
Joining, by = "factor"
Joining, by = "factor"

efa13_loadings <- efa13$loadings[] %>%
  data.frame() %>%
  rownames_to_column("question") %>%
  gather(factor, loading, -question) %>%
  mutate(scale = gsub("p7_", "", question),
         scale = gsub("_.*$", "", scale),
         factor = factor(factor, levels = paste0("MR", 1:13)))

efa13_loadings_dom <- efa13_loadings %>%
  group_by(scale, question) %>%
  top_n(1, abs(loading)) %>%
  ungroup() %>%
  arrange(factor, desc(abs(loading)))

efa13_loadings_dom %>%
  count(factor, scale) %>%
  complete(factor, scale, fill = list(n = ".")) %>%
  spread(scale, n) %>%
  kable() %>%
  kable_styling()

factor	abs	dse	exsen	hthk	mm	por	se	unev	wob
MR1	.	13	.	.	.	3	11	.	.
MR2	32	.	.	.	.	.	.	.	.
MR3	2	.	.	.	.	.	.	.	6
MR4	.	.	.	.	14	3	.	.	.
MR5	.	.	.	.	.	.	11	.	.
MR6	.	.	.	.	.	.	.	.	6
MR7	.	1	.	9	.	.	.	.	.
MR8	.	.	.	.	.	.	.	6	.
MR9	.	.	.	9	.	.	.	.	.
MR10	.	.	.	.	3	.	.	.	.
MR11	.	.	.	.	.	10	.	.	.
MR12	.	.	7	.	.	.	.	.	.
MR13	.	.	.	.	.	.	.	.	.

Quick interpretations:

• MR1 = spiritual experience I (dse/se)
• MR2 = absorption
• MR3 = other personality (wob/abs)
• MR4 = martha mary
• MR5 = spiritual experience II (se)
• MR6 = sense of control (wob)
• MR7 = need for cognition I (hthk)
• MR8 = launay-slade (unev)
• MR9 = need for cognition II (hthk)
• MR10 = other porosity
• MR11 = porosity
• MR12 = sheep-goat (exsen)
• [MR13 = NULL]

Focus on items that from DSE, SE, Porosity, & MM scales:

efa13_loadings_dom %>% 
  mutate(dom = "bold") %>% 
  select(-loading) %>%
  full_join(efa13_loadings) %>%
  full_join(efa13_loadings_dom %>%
              arrange(desc(factor), abs(loading)) %>%
              mutate(order = 1:nrow(.)) %>%
              select(scale, question, order)) %>%
  mutate(dom = ifelse(is.na(dom), "plain", dom)) %>%
  filter(scale %in% c("dse", "se", "por", "mm")) %>%
  ggplot(aes(x = factor, y = reorder(question, order),
             fill = loading, label = format(round(loading, 2), nsmall =))) +
  facet_grid(scale ~ ., scales = "free", space = "free") +
  geom_tile(color = "black") +
  geom_text(aes(size = dom, fontface = dom)) + 
  scale_fill_distiller("Factor loading", palette = "RdYlBu", limits = c(-1, 1),
                       guide = guide_colorbar(barheight = 20, barwidth = 1)) +
  scale_x_discrete(position = "top") +
  scale_size_manual("Dominant factor?", 
                    values = c(3, 2), labels = c("dominant", "not dominant")) +
  theme_minimal() +
  theme(axis.title = element_blank())

Joining, by = c("question", "factor", "scale")
Joining, by = c("question", "scale")

Minimizing BIC

VSS(d1)

Very Simple Structure
Call: vss(x = x, n = n, rotate = rotate, diagonal = diagonal, fm = fm, 
    n.obs = n.obs, plot = plot, title = title, use = use, cor = cor)
VSS complexity 1 achieves a maximimum of 0.72  with  1  factors
VSS complexity 2 achieves a maximimum of 0.8  with  3  factors

The Velicer MAP achieves a minimum of 0  with  8  factors 
BIC achieves a minimum of  -42163.96  with  6  factors
Sample Size adjusted BIC achieves a minimum of  -12097.56  with  8  factors

Statistics by number of factors 

Minimizing BIC suggests retaining 6 factors.

efa6 <- fa(d1, 6, rotate = "varimax")

heatmap_fun(efa6, factor_names = paste0("MR", 1:6))

the condition has length > 1 and only the first element will be usedJoining, by = "capacity"
Joining, by = "factor"
Joining, by = "factor"

efa6_loadings <- efa6$loadings[] %>%
  data.frame() %>%
  rownames_to_column("question") %>%
  gather(factor, loading, -question) %>%
  mutate(scale = gsub("p7_", "", question),
         scale = gsub("_.*$", "", scale),
         factor = factor(factor, levels = paste0("MR", 1:6)))

efa6_loadings_dom <- efa6_loadings %>%
  group_by(scale, question) %>%
  top_n(1, abs(loading)) %>%
  ungroup() %>%
  arrange(factor, desc(abs(loading)))

efa6_loadings_dom %>%
  count(factor, scale) %>%
  complete(factor, scale, fill = list(n = ".")) %>%
  spread(scale, n) %>%
  kable() %>%
  kable_styling()

factor	abs	dse	exsen	hthk	mm	por	se	unev	wob
MR1	.	14	.	.	.	1	10	.	.
MR2	2	.	.	6	.	.	.	.	6
MR3	32	.	4	4	.	.	.	1	.
MR4	.	.	3	.	17	15	.	.	.
MR5	.	.	.	.	.	.	12	5	.
MR6	.	.	.	8	.	.	.	.	6

Quick interpretations:

• MR1 = spiritual experience I (dse/se)
• MR2 = other personality (hthk/wob/abs)
• MR3 = absorption (plus some exsen/hthk/unev)
• MR4 = porosity (mm/por)
• MR5 = spiritual experience II (se/unev)
• MR6 = control (hthk/wob)

Focus on items that from DSE, SE, Porosity, & MM scales:

efa6_loadings_dom %>% 
  mutate(dom = "bold") %>% 
  select(-loading) %>%
  full_join(efa6_loadings) %>%
  full_join(efa6_loadings_dom %>%
              arrange(desc(factor), abs(loading)) %>%
              mutate(order = 1:nrow(.)) %>%
              select(scale, question, order)) %>%
  mutate(dom = ifelse(is.na(dom), "plain", dom)) %>%
  filter(scale %in% c("dse", "se", "por", "mm")) %>%
  ggplot(aes(x = factor, y = reorder(question, order),
             fill = loading, label = format(round(loading, 2), nsmall =))) +
  facet_grid(scale ~ ., scales = "free", space = "free") +
  geom_tile(color = "black") +
  geom_text(aes(size = dom, fontface = dom)) + 
  scale_fill_distiller("Factor loading", palette = "RdYlBu", limits = c(-1, 1),
                       guide = guide_colorbar(barheight = 20, barwidth = 1)) +
  scale_x_discrete(position = "top") +
  scale_size_manual("Dominant factor?", 
                    values = c(3, 2), labels = c("dominant", "not dominant")) +
  theme_minimal() +
  theme(axis.title = element_blank())

Joining, by = c("question", "factor", "scale")
Joining, by = c("question", "scale")

efa6_loadings_dom %>% 
  mutate(dom = "bold") %>% 
  select(-loading) %>%
  full_join(efa6_loadings) %>%
  full_join(efa6_loadings_dom %>%
              arrange(desc(factor), abs(loading)) %>%
              mutate(order = 1:nrow(.)) %>%
              select(scale, question, order)) %>%
  mutate(dom = ifelse(is.na(dom), "plain", dom)) %>%
  # filter(scale %in% c("dse", "se", "por", "mm")) %>%
  ggplot(aes(x = factor, y = reorder(question, order),
             fill = loading, label = format(round(loading, 2), nsmall =))) +
  facet_grid(scale ~ ., scales = "free", space = "free") +
  geom_tile(color = "black") +
  geom_text(aes(size = dom, fontface = dom)) + 
  scale_fill_distiller("Factor loading", palette = "RdYlBu", limits = c(-1, 1),
                       guide = guide_colorbar(barheight = 20, barwidth = 1)) +
  scale_x_discrete(position = "top") +
  scale_size_manual("Dominant factor?", 
                    values = c(3, 2), labels = c("dominant", "not dominant")) +
  theme_minimal() +
  theme(axis.title = element_blank())

Joining, by = c("question", "factor", "scale")
Joining, by = c("question", "scale")

Weisman et al. (2017) retention criteria

reten_fun(d1, "none")

[1] 2

The retention criteria employed in Weisman et al. (2017) suggest retaining 2 factors.

efa2 <- fa(d1, 2, rotate = "varimax")

heatmap_fun(efa2, factor_names = paste0("MR", 1:2))

the condition has length > 1 and only the first element will be usedJoining, by = "capacity"
Joining, by = "factor"
Joining, by = "factor"

efa2_loadings <- efa2$loadings[] %>%
  data.frame() %>%
  rownames_to_column("question") %>%
  gather(factor, loading, -question) %>%
  mutate(scale = gsub("p7_", "", question),
         scale = gsub("_.*$", "", scale),
         factor = factor(factor, levels = paste0("MR", 1:2)))

efa2_loadings_dom <- efa2_loadings %>%
  group_by(scale, question) %>%
  top_n(1, abs(loading)) %>%
  ungroup() %>%
  arrange(factor, desc(abs(loading)))

efa2_loadings_dom %>%
  count(factor, scale) %>%
  complete(factor, scale, fill = list(n = ".")) %>%
  spread(scale, n) %>%
  kable() %>%
  kable_styling()

factor	abs	dse	exsen	hthk	mm	por	se	unev	wob
MR1	2	13	4	8	17	16	19	.	5
MR2	32	1	3	10	.	.	3	6	7

Quick interpretations:

• MR1 = porosity + spiritual experience? (“spiritual”?)
• MR2 = absorption + personality? (“secular”?)

Focus on items that from DSE, SE, Porosity, & MM scales:

efa2_loadings_dom %>% 
  mutate(dom = "bold") %>% 
  select(-loading) %>%
  full_join(efa2_loadings) %>%
  full_join(efa2_loadings_dom %>%
              arrange(desc(factor), abs(loading)) %>%
              mutate(order = 1:nrow(.)) %>%
              select(scale, question, order)) %>%
  mutate(dom = ifelse(is.na(dom), "plain", dom)) %>%
  filter(scale %in% c("dse", "se", "por", "mm")) %>%
  ggplot(aes(x = factor, y = reorder(question, order),
             fill = loading, label = format(round(loading, 2), nsmall =))) +
  facet_grid(scale ~ ., scales = "free", space = "free") +
  geom_tile(color = "black") +
  geom_text(aes(size = dom, fontface = dom)) + 
  scale_fill_distiller("Factor loading", palette = "RdYlBu", limits = c(-1, 1),
                       guide = guide_colorbar(barheight = 20, barwidth = 1)) +
  scale_x_discrete(position = "top") +
  scale_size_manual("Dominant factor?", 
                    values = c(3, 2), labels = c("dominant", "not dominant")) +
  theme_minimal() +
  theme(axis.title = element_blank())

Joining, by = c("question", "factor", "scale")
Joining, by = c("question", "scale")

Four-factor solution

Just for the hell of it, here’s a four-factor solution.

efa4 <- fa(d1, 4, rotate = "varimax")

heatmap_fun(efa4, factor_names = paste0("MR", 1:4))

the condition has length > 1 and only the first element will be usedJoining, by = "capacity"
Joining, by = "factor"
Joining, by = "factor"

efa4_loadings <- efa4$loadings[] %>%
  data.frame() %>%
  rownames_to_column("question") %>%
  gather(factor, loading, -question) %>%
  mutate(scale = gsub("p7_", "", question),
         scale = gsub("_.*$", "", scale),
         factor = factor(factor, levels = paste0("MR", 1:4)))

efa4_loadings_dom <- efa4_loadings %>%
  group_by(scale, question) %>%
  top_n(1, abs(loading)) %>%
  ungroup() %>%
  arrange(factor, desc(abs(loading)))

efa4_loadings_dom %>%
  count(factor, scale) %>%
  complete(factor, scale, fill = list(n = ".")) %>%
  spread(scale, n) %>%
  kable() %>%
  kable_styling()

factor	abs	dse	exsen	hthk	mm	por	se	unev	wob
MR1	.	14	.	1	.	.	19	.	.
MR2	33	.	6	4	.	.	3	6	.
MR3	1	.	.	7	.	.	.	.	9
MR4	.	.	1	6	17	16	.	.	3

Quick interpretations:

• MR1 = spiritual experience
• MR2 = absorption (+ secular experience?)
• MR3 = controls
• MR4 = porosity (+ personality?)

Focus on items that from DSE, SE, Porosity, & MM scales:

efa4_loadings_dom %>% 
  mutate(dom = "bold") %>% 
  select(-loading) %>%
  full_join(efa4_loadings) %>%
  full_join(efa4_loadings_dom %>%
              arrange(desc(factor), abs(loading)) %>%
              mutate(order = 1:nrow(.)) %>%
              select(scale, question, order)) %>%
  mutate(dom = ifelse(is.na(dom), "plain", dom)) %>%
  filter(scale %in% c("dse", "se", "por", "mm")) %>%
  ggplot(aes(x = factor, y = reorder(question, order),
             fill = loading, label = format(round(loading, 3), nsmall =))) +
  facet_grid(scale ~ ., scales = "free", space = "free") +
  geom_tile(color = "black") +
  geom_text(aes(size = dom, fontface = dom)) + 
  scale_fill_distiller("Factor loading", palette = "RdYlBu", limits = c(-1, 1),
                       guide = guide_colorbar(barheight = 30, barwidth = 1)) +
  scale_x_discrete(position = "top") +
  scale_size_manual("Dominant factor?", 
                    values = c(3, 2), labels = c("dominant", "not dominant")) +
  theme_minimal() +
  theme(axis.title = element_blank())

Joining, by = c("question", "factor", "scale")
Joining, by = c("question", "scale")

# EXTRA
# hierarchical clustering
clust <- d1 %>% t() %>% dist() %>% hclust()
clust %>%
  ggdendrogram(rotate = T) +
  theme_minimal() +
  theme(panel.grid.major.y = element_blank(),
        panel.grid.minor.y = element_blank()) +
  labs(title = "Hierarchical agglomerative clustering", 
       subtitle = "Complete linkage (default for stats::hclust() function)",
       y = "Height", x = "Question")